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Abstract 

We present a non-standard eigenvalue problem that arises in the linear stability of a three-layer 
Hele-Shaw model of enhanced oil recovery. A nonlinear transformation is introduced which allows refor¬ 
mulation of the non-standard eigenvalue problem as a boundary value problem for Kummer’s equation 
when the viscous profile of the middle layer is linear. Using the existing body of works on Kummer’s 
equation, we construct an exact solution of the eigenvalue problem and provide the dispersion relation 
implicitly through the existence criterion for the non-trivial solution. We also discuss the convergence 
of the series solution. It is shown that this solution reduces to the physically relevant solutions in two 
asymptotic limits: (i) when the linear viscous profile approaches a constant viscous profile; or (ii) when 
the length of the middle layer approaches zero. 
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1 Introduction 


The flow of two immiscible fluids through porous media arises in many important industrial and natural 
situations such as secondary oil recovery, ground water remediation, and geological CO 2 storage. Such flows 
are known to be potentially unstable, especially when the displacing fluid is more viscous than the displaced 
one. There exist some similarities between porous media and Hele-Shaw flows (i.e. flow in a Hele-Shaw cell, 
see below); for example the pressure drop in both such flows are governed by Darcy’s law for single fluid flow. 
Due to this and the fact that it is significantly easier to study Hele-Shaw flows theoretically, numerically, 
and experimentally, there have been numerous theoretical and numerical studies even for Hele-Shaw flow 
of two immiscible fluids since the early 1950s, starting with the work of Saffman and Taylor [10]. There 
are many review articles on such studies, for example see BM- These studies were originally motivated by 
displacement processes arising in secondary oil recovery, even though these studies have much wider appeal 
in the sciences and engineering. In the late 1970s, tertiary displacement processes involved in chemical 
enhanced oil recovery generated interest in three-layer and multi-layer Hele-Shaw flows (see [U [S] [5] [7] ). 

In this paper, we first briefly derive the non-standard eigenvalue problem. This eigenvalue problem has 
been derived earlier by the first author and his collaborators; for example see |5]. But the difference is that 
the derivation presented here is more general and shows how to generate higher order correction terms if 
necessary in order to study the effect of nonlinear terms that may dominate the dynamics, particularly in 
view of the sensitivity of fingering problems to finite amplitude perturbations. However, we do not study 
or discuss such nonlinear effects in this paper which will be taken up in the future as it falls outside the 
scope of this paper. We then analytically study this non-standard eigenvalue problem using non-linear 
transformation for the case when the viscous profile of the middle layer is linear. We will see below that this 
case is relatively hard to study in comparison to the case when the viscous profile is exponential which we 
have recently addressed in [6]. 



Figure 1: Three-layer rectilinear Hele-Shaw flow in which the middle layer has a smooth viscous profile. The 
physical set-up as well as the smooth viscous profile of the middle layer are shown in this figure. 

The physical set-up consists of rectilinear motion of three immiscible fluids in a Hele-Shaw cell which is 
a device separating two parallel plates by a distance b (see Fig. [T|). The fluid in the extreme left layer Ri 
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with viscosity /xi extends up to a; = —oo, the fluid in the extreme right layer i ?2 with viscosity ^2 > Mi 

extends up to a; = 00 , and the fluid in the middle-layer Rj of finite length L has a smooth viscous profile 

with viscosity increasing in the direction of displacement. The interfacial tensions of the leading and the 
trailing interfaces are given by T and S respectively. It is well established that this Hele-Shaw flow is similar 
to flow in homogeneous porous media with equivalent permeability 6^/12. Without any loss of generality, 
we take this to be one below. The mathematical model considered here consists of conservation of mass, 
Darcy’s law and advection equation for viscosity. Thus we have 

V-u = 0, V(a;,y)eR^ (1) 

{ Mi; if {x,y) € i?i 

M(x,?/,t); if (x,?/) G i?/ (2) 

M 2 ; if {x,y) G i ?2 

y,t -\- uyix + vy,y =0; \f x € Rj and y G K. (3) 

Due to the continuity equation, we can define the stream function ip = tpix,y,t) such that u = ipy and 
V = —ipx- This then implies that 

Px = -pipy, Py = pipx, and Mt + IpyPx - ipxPy = 0. (4) 

Since u = {Uq, 0) when -I- —>■ oo, we consider a small perturbation of the basic scalar fields ipo,PQ and 

Mo of the form 

Ip ='ip{x,y,t) = Uoy + eip{x,y,t) 

P =p{x,y,t)=po{x,t)+ep{x,y,t) > 

M = p{x,y,t) = po{x,t) + ep{x,y,t) 

Substituting into the original equations, we get the following 0{e°) and 0{e^) equations. 

0(e°) equations: 

POx = —UqPo 

POy = 0 > 

Pot + Uopox = 0 

/ 

These equations provide the basic solution given by 

po = po{x,t) = po{x - Uot) I 

Pq = PQ{x,t) =-Uo po{s-Uot)ds > 

Xq 

u=(C/o,0) 

where po{x — Uot) is an arbitrary function of ^ = x — Uot, meaning the viscous profile is fixed with respect 
to a moving frame moving at a constant velocity (Uq, 0). 

O(e^) equations: 

^ ' 

Px = -UoP - Polpy 

Py = Polpx > 

Pt + UoPx + POxIpy = 0. 
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Now, introducing the moving frame change of variables, namely ^ = x — U^t, y 
following system of equations. 

^ ' 

Py = 

Mt + = 0- 


y, t = t, we get the 


(5) 


Taking cross derivatives of the hrst two equations with respect to y and ^ respectively and then subtracting 
the resulting equations from each other gives + U^jly + = 0. This combined with the equation 

(inilg leads to 

- Uoyo^yy + Po(A^)t = 0. 


Using the ansatz ■0 = in the above equation together with the appropriate boundary conditions 

(see i) give the following eigenvalue problem. 


fJ-oiOihi - + miOfi + = 0 

= fi-L) jpiA: + —^[mi - P(l"(-^)] + — j 

Mo (0)/c(0) = /(O) |-M 2 fc+ ^^[M 2 -M0(O)] - 


where the viscous prohle of the middle layer, namely, yo = ^o(C) is an arbitrary function. This is a non¬ 
standard eigenvalue problem in that the spectral parameter l/cr appe ars in the equation as well as in 
the boundary conditions. Recently, this problem has been numerically solved by Daripa [3] for a constant 
viscous profile and by Daripa & Ding [4] for non-constant viscous profiles /ro(^) to determine the most 
optimal profile, i.e., the least unstable profile. This problem has been too difficult to solve analytically for 
non-constant prohles. Progress made in this direction for the linear viscous profile is presented below. 

In this paper we consider a linear viscous profile for the intermediate fluid region given by 


Mo(C) = + /3 for - L < ^ < 0 


where 

(M2 - Ml) - ('^1 + >^2) Mo( 0 ) -/ro(-T) ^ ^ ^ 

a = - - -= - - -, 0 = ^2 - M2 = Mo( 0 ), 

and Ji = fj,Q{—L) — yi, J 2 = tJ -2 — Mo(0) are jump discontinuity values at the interfaces ^ = —L and ^ = 0, 
respectively. In the left region the problem reduces to = 0,lim^_>_oo /(■?) = O 7 which has solution 

/(^) = /(—for ^ < —L. In the right region the problem reduces to — k'^f = 0, hm^_>oo /(C) = 0; 
which has solution /(C) = /(0)e“^^ for C > 0. In the intermediate region the problem reduces to 


[ (aC + /3)/cc+a/4 + fc2(aA-(aC + /3))/ 
< no{-L)f^{-L) 

[ Mo(0)/a0) 


0, -L < C < 0 
(ai(fc)A + a2(fc))/(-L) 
(/3i(fc)A-/?2(fc))/(0) 


where A = ^ is the spectral parameter and 


ai(fc) = + ^^(Mi - Mo(-i)), 

/3i(fc) = -^ + ^'(M2-Mo(0)), 


( 6 ) 
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a2{k) = y,ik 
/32(fc) = yi2k 


(7) 








2 Solution to the Eigenvalue Problem (E]) via Kummer’s Equation 


We introduce the nonlinear transformation and change of variables given by 


fiO = w = —2k ^ < 0. 


( 8 ) 


After some manipulation of the eigenvalue problem (|6]) using the above transformation, we obtain the 
following eigenvalue problem for the Kummer’s equation 


wZww + {b — w)zyj — az = 0, W 2 < w < wi 
r]iz{wi) + (l)iz'{wi) = 0, 

r]2Z{w2) + (I>2z'iw2) = 0, 

where a prime denotes derivative, 

& = 1, a = -(1 + fcA), W 2 = w{^ = 0) = —2k \ wi = w(^ = —L) = W 2 + 2kL < 0, 

2 a 

r]i = Sk^/a, (j)i=2k^i^ r ]2 = 2^.2k + Tk'^/a, and 4>2 = —2kiJ,2- 


(9) 


( 10 ) 


The eigenvalue problem Q is a regular two point boundary value problem for each wave number k. One 
solution zi of the Kummer’s equation is given by 

9 rj 

aiw a2W^ asW^ anW^ 

Zi[a, 1, w) — flo T /, INo H——I— /r.ixr, + ■ H——rr;^—h ■ 


(1!)2 (2!)2 (3!)2 


( n !)2 


where oq = 1 and a„ = a(a + l)(a + 2) ■ ■ ■ (a + n — 1) for n = 1 , 2 ,.... 

This is an analytic solution. It is easily seen that the derivative of this solution which we will need below 
for the dispersion relation is given by 


, , , , ai 2 a 2 W 3a3W^ nanw’^ ^ (n + l)an+iw'’ 

Zi(a, l,w) =H-^- . + ■ 


+ • • • = azi(a + 1, 2, w). 


(1!)2 (2!)2 (3!)2 (n!)2 (n+l)!2 

The linearly independent second solution is easily constructed by the method of Frobenius. Avoiding all 
the details, the second solution Z 2 is given by 


Z2(a, 1, w) = ^ 1 + 


ai/2W 


+ 


0-3/21 


O 3 / 2 I 


+ ••• + 


+ ■ 


1-3-1! 1-3-5-2! 1-3-5 •7-3! 1 • 3 • 5 • 7 • • • (2n + 1) ■ n! 

where = (a + • • • (a + |)(a + |)(a + A), n = 1,2,3,.... Its derivative which we will 
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need below is then given by 


^2 


/ 1 \ 1 ^ 1 \ , 1/2 j 2 ai /2 , 2 • 2^a3/2W 

{a,l,w) =—Z 2 {a,l,w) + w ' -1--h 


n2"a. 


3! 


5! 


7! 


+ ••• + 


(2n + l)! 


+ ■ 


The general solution of the Kummer’s equation is then given by 

z{w) = Cl2:1(0, 1, w) + 022:2(0, 1, w), 


where ci and C 2 are arbitrary constants. Substituting the general solution into the two boundary conditions 
of the eigenvalue problem ([9]), we obtain the following linear system of equations for ci and C 2 . 

{ 77 lZl(wi) + 4 >iz[{wi)]ci + {7712:2(101) + 0iZ2(wi)}c2 = 0 , 

{■q2Zl{w2) + (j)2z'i{w2)}ci + {t]2Z2{w2) + «!)24(^2)}c2 = 0. 
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Therefore, for a non-trivial solution we have 


rjiziiwi) + <piz[{wi) 77122 (^ 1 ) + 'piZ 2 {wi) 
r]2Zi(w2) + 4'2 z[{w 2) r]2Z2iw2) + 4'2Z2{w2)- 


= 0 . 


( 11 ) 


This formally gives the dispersion relation (T{k) in terms of the problem data: S, T, /ii, /i 2 , T and Uo- 

In terms of the original variables / and the fundamental solutions fi{^) and / 2 (C) are then given by 
(see (HI)) 


. fef r. «i2fc(e + f ) , a2{2kr{^ + as{2kn^ + f )3 ^ 

M0=e |1--^yy^-+ 


+ (-l)^ 


Ii2 


/.({)=='•* {-2*, (£ + ^)}‘'"|l 


(n!) 

Ai2\'l^/^/, ai/22fc(.^ + ^) a3/2(2fc)^(^ + 


( 12 ) 


1 ■ 3 ■ 5 • 7 - (3!) 


+ ••• + (- 1 )^ 


l-3-(l!) l-3-5-(2!) 

a.„_.(2fc)"(e + ^)^ 


l-3-5-7---(2n + l)(n!) 


(13) 


where 


ao = 1 

a„ = a(a + l)(a + 2) • • • (a + n — 1) 

a 2 n-i = (a + l/2)(a + 3/2)(a + 5/2) • • • (a + {a + , 1, 2, 3,4, 5,... 

a_i = 1 
2 

Also, recall that ^ = x — Uot, a = ( 7^2 — Mi)/!^ and a = {1 + kUo/cr)/2. Noticing that both series are centered 
at ^ = — (/r 2 /a) and applying the ratio test for series for /i(^) above, we have 


((n + 1)!)2 


(n!)2 


(a + n)(2fc)(C + ^) 


(n + l)2 


(a + n){2k) 


(n + l)2 






0 as n —?► 00 , 


for any ^ fixed. Thus the series for /i(^) converges absolutely V ^ — (/r 2 /a), but the series evaluated at 

^ = —(^ 2 / 0 ) reduces to 1. This implies that the radius of this series is 00 . Hence 

(-l)"a„(2A:)’"(^+ ^)’" 
lim 1— ^ \ ^ = 0, 

n->oo (n!) 

and since it is an alternating singular series, the error En^ in approximating fi by terms up to fc" is smaller 
than the last neglected term, namely 

a„+i(2A:)-+i(^+i^)-+i 




((n + l)!)2 

for a fixed Similarly, applying the ratio test to the series for /2(0; 


(2n + 3)! 


I M2'\n 
'a ■' 


(2n+ 1)! 


2(a+2i^)(2fc)(^+i^) 


(2ti T 2) (2/1 T 3) 


0, as n —?> 00 , 
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for any ^ fixed. Thus the series inside the brackets converges absolutely ^ — (/X 2 /Q;), but the series inside 
the brackets evaluated at ^ alpha reduces to 1. Therefore, the radius of convergence of the series 

within the brackets is 00 . Hence 

lim -;-= 0. 

n —>00 (2n + 1)! 


Notice that f2{0 has a branch point at ^ In any case, since the series within the brackets is an 

( 2 ) 

alternating sign series, if we truncate it, the error En is smaller than the last neglected term, i.e.. 






(2n + 3)! 


The general solution of the ODE ([6])^ is then given by /(^) = Ci/i(^) + C2f2{0- The boundary values of 
/(^) follow from (|T^ and (fT^ whi ch are now given by 


n—0 

CX) 


(niy 


/l( 0 ) = ^ ^ ^ 

n—0 


{n\y 


/2(-L) = e-'=^|-2fcT 

/2(0) = |-2fcL 


Ml 


M2 - Ml 


1 / 2 ^ (-l)"(22A:Lra,„_,(^r 


E 

n—O 


(2n + l)! 


M2 


M2 - Ml 


1 / 2 ^ (_l)n(22/cL)"a,„_, 


E 

n—0 


(2n+l)! 


Substituting these in the boundary conditions ([51)3, we obtain the following system of equations for the 
constants Ci and C 2 . 


— fiif[{—L)}ci + {Aif 2 {—L) — pif 2 {—L)}c 2 — 0, 
{^ 2 / 1 ( 0 ) - M2/l(0)}ci + {^ 2 / 2 ( 0 ) - /i2/^(0)}c2 = 0, 


where Ai = {pik + and A 2 = {—p 2 k + 2-^}. For the existence of nontrivial solutions, we then have 


A,h{-L) - fiifii-L) A,f2i-L) - pif^{-L) 
212/1(0) - P 2 fm ^2/2(0) - P 2 fm 


which gives us the dispersion relation in the form: 9{<t, k) = 0. Because of the nature of the series solutions 
given above, it is not possible to give this dispersion relation explicitly. 


3 Limiting Cases 

There are an infinite number of eigenvalues a (recall A = U/cr) which can be ordered: CTmax = cti > 172 > 

. > (Too —5" 0. We know that these infinite number of eigenvalues should reduce to (i) only two in the 

limit q: —?> 0 corresponding to the constant viscosity of the intermediate layer fluid (see Daripa [3]); (ii) only 
one in the limit L —>■ 0 (see Saffman & Taylor [TU], Daripa [2]) and (iii) only two in the limit of L —>■ 00 
(see Daripa HD- In fact, we also know the eigenvalues in these limiting cases from the pure Saffman-Taylor 


7 



















growth rate of individual interfaces. These results by no means are transparent from the solutions of the 
eigenvalue problem (jS]) given in the previous section. Below, we show how to recover these limit solutions 
(eigenvalues) from the infinite number of eigenvalues for the linear viscous profile. 


3.1 Constant viscosity case: a = 0. 

In this case, the eigenvalue problem ([S]) reduces to 

[ = 0, -T<e<0, 

I ^o(-L)/ 5(-L) = (ai(fc)A + a2(A:))/(-L), (14) 

[ Mo(0)/c(0) = (^i(fc)A-/32(A:))/(0), 

In this case, the change of variable introduced previously, namely w = —2k , which converts the 

equation (|6|)j^ to Rummer’s equation, i s not well-defined. Therefore we work with the boundary value 
problem (fHl) . Now, consider the general solution of the ODE in (imi 

f{0 = Ah{0+Bf2{0. 

such that 

r /i(-L) = i, /{(-i) = o, 

\ h{-L) = 0, /^(-L) = 1. 

Therefore, we get 

/i(0 = cosh(fc(^ -k L)), and /2(0 = ^ 

We search for a solution of the boundary value problem (imi of the form 

/(^; A) = Alcosh(fc(e + £)) + g smh(fc(g + L)) 

k 

To find a solution of (IT4l) of this form, we start determining the coefficients using the shooting technique 
such that the boundary condition at ^ = —L is satisfied. Obviously, the coefficients A and B depend on the 
parameter A. Then, we find A in such a way that the solution satisfies the boundary condition at = 0. 
Hence, we look for A and B such that 


f{-L;X)=fj,o{-L), and /^(-L; A) = ai(fc)A-b a 2 (^)- 
Then it follows directly from (fTSl) that A = iJ,o{—L) and B = ai{k)X -b a 2 {k). Therefore 

fiC, X) = Mo(-i) cosh(fc(e + L)) + (ai(fc)A + 


(17) 


satisfies the ODE in (ED and the boundary condition at ^ = —L. Erom these it follows that the spectrum 
of problem ED can be studied using the following algebraic equation (see EDa) 


Mo(0)4w^ = /5i(fe)A-/32(fc), 


(18) 


/(0;A) 

where f{^]X) is the function defined in (fT71) . Evaluating /(0;A) and f^{0-X) from (fT71) and substituting 
directly in (fTsP one obtains 


fc^o(O) < 


f fio{—L) sinh(fcL) -b (ai(A:)A -b 0 : 2 (fc)) *^^^ ^ ^ 


I cosh(fcL) -b {ai{k)X + 0 ^ 2 (fc)) 

k 


= A(fc)A-/32(fc). 


(19) 










Then, taking i —> O'*" one obtains 


which is equivalent to the equation 


{fio{0)ai{k) - ^o{-L)l3i{k)) A = - {fio{0)a2{k) + fio{-L)l32{k)). 
Since HoiO = constant {fi or /3), it follows that 


(ai(fc) - /?i(fc))A = -{a2{k) + (32{k)). 


Now, using the definition of the coefficient ai(fe), a 2 {k), /3i(fe), I32{k) given in ([7]) we have 


Tk^ \ 





A = -k{fj.i + fi2) 


from which it follows that 

^ _ Uok{n2 - A^i) - k^{S + T) 

(/^i +M2) 

which is the formula for the growth rate of an interface with surface tension (S' + T), which is what should 
be expected in this limit. Thus we recover the classical formula for the growth rate in this limit. 

To take the limit when L —>■ 00 , we go back to equation m and write it as follows 


kfioi.0) 


^o{—L) tanh(/cT) + (ai(fc)A + a 2 (k))- 


Ho{-L) + (Q;i(fc)A + a2(fc)) 


tanh(fcT) 


/?i(fc)A-/32(fc). 


Now, taking the limit when L 00 , we obtain fc^o(O) = /3i(/c)A — (32(k). Using a = Uq/X and expressions 
for the coefficients from (l7|)2, we obtain 


( -^ + kHf^2-M0)) \ 

k{^i2 - ^lo{0)) )' 


Finally, since fio{Q) = constant or /3) it follows that 


Tk^ 

a =-h Uok 

M2 + M 


/ M2 - 
VM + M2/ 


which gives the classical formula for Saffman-Taylor instability of the leading interface. Similarly, we can 
recover the the classical formula for Saffman-Taylor instability of the trailin g interface by reversing the 
shooting technique (see after dTSll b i.e., first find the solution which is analogous to (ITTl) but satisfies the 
boundary condition at ^ = 0 instead and then shoot to satisfy the boundary condition at ^ = —L (i.e., 
replace (1181) by a similar formula derived from the boundary condition at ^ = — L and follow the procedure). 


3.2 Linear viscosity case: a > 0. 

In this section, we study asymptotic limits (L —>■ 0 and L 00 ) of the solutions to the eigenvalue problem (O. 
To this end, we consider the following form of two linearly independent solutions of Rummer’s equation @ 1 . 
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These are convenient for the asymptotic analysis presented below. 


e^U{l-a,l,-w) = - 


r(l-a) 


M(1 — a, 1, —w) ln(—ii;)+ 


VJ /1 

- WTTTT E -77^ -«) + *)- 2'P(1 + *)) {-wT 


(20) 


r(l-a)^ (*!)2 

where 'I'(s) is Euler’s digamma function (See Abramowitz [1], Chapter 13). 

To this end, we follow the steps presented in the previous section [XT] for the particular case = 

constant (/i or 0). From the transformation in ([8]), it follows that 

f{0A) = {Azi{w) + Bz2{w))e^^ ( 21 ) 

where 


is the general solution of the ODE 

zi{w) = C\M{a,l,w) + Die^U{l — a0,—w) 
Z 2 (w) = C 2 M{a,l,w) + D 2 e'^U{l — a,l,—w) 

and Cl, Di, C 2 , D 2 are chosen such that 


Zi{wi) = 1, z[{wi) = 0 1 

Z2iwi)=0, Z2{wi) = l J 

where a = (1 + kX)/2. Substituting (|22|) in the boundary conditions 
systems of equations 

( M(a, 1, rci) e“"^t/(l — a, 1, —wi) \ f Ci 


( 22 ) 


(23) 


M'{a,l,wi) (e"'f/(l — a, 1, —w)(„ 


.Di, 


we obtain the following linear 

(24) 


/ M(a, 1, wi) e’"it/(l - a, 1, -wi) | \ / 0\ 

l^M'(a,l,a;i) (e“t/(l - a, 1,-H)^/ W ~ VV 

Solving the above two systems and using the relations (see Abramowitz [1], Chapter 13) 

M'{a,l,w) = aM(a + 1,1 + 1, ic), 1 

U'{l — a,l,—w) = —(1 — a)[/(l + (1 — a), 1 + 1, —w)(—1). J 


(25) 


(26) 


we obtain 


Cl = e“^(c/(^E,i,-«;i)+l-Ec(^E,2,-«;i) )/lE{l,2} 


1 - kX 


C 2 = 


D 2 = M( i^,l,u;i )/lE{l,2} 


,l,-u;i /1T{1,2} 


(27) 


where 1E{1, 2} is the determinant of the coefficient matrix of the system 

Similar to the procedure of the previous section [3T] we find A and B so that f{—L;X) = fio{—L) and 
f'i—L; X) = ai{k)X + a 2 {k). Therefore, it follows from (1^ and (051) that 

Ae~'"^ = ^o(-T), 

Ake~^^ —2kBe~^^ = ai{k)X + a2{k), 
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and therefore A = fj,o{—L)e’^^ and B = — ^ '^dfc)-^+a 2 (fc)-fcAto(-^) j ^kL_ Substituting these constants in the 
function /(^; A) defined by (I^Tl) . we obtain a solution of the ODE that satisfies the boundary condition at 
^ = —L of the eigenvalue problem (O. Since A and B depend on the spectral parameter A, it follows that 
the eigenvalues of the problem ([S]) can be obtained by studying the following algebraic equation which is a 
reformulation of the boundary condition at ^ = 0 of the eigenvalue problem (|6l) . 

MO)^0^ = Mk)X-P2ik). ( 28 ) 

Since the right-hand side of the above equation does not depend on D, we need to study the asymptotic 
limits (L —?> 0 and L —>■ oo) of the lefthand side of (1^51) . Notice that the expression /t(0; A)//(0; A) above is 
given by (see (ED) 

/g(0;-^) = r, /^i _ + Bz!^{w 2 ) \ , . 

/(0;A) Iv Az,{w2) + Bz2iw2)) ' ^ ’ 

Therefore, we first find the asymptotic approximations for 21 (^ 2 ), ^ 2 (^ 2 ), z((w 2 ), and - 22 (^ 2 ) in both cases 
below before estimating the ratio /^(O; A)//(0; A) using (E^D for its use in (ED- Below, we write wi = diL 
and W 2 = d 2 L where di and d 2 are given by (see (HHD), 

=-2kfj.o{-L)/{fio{0) - d-o{-L)), and d 2 =-2k{fio{0))/{fJ-o{0) - do{-L)). (30) 


First case (When L —>■ 00 ): It follows from Abramowitz and Stegun [I] that 

M{a,l,djL) = ^^(-djL)"“(l-H0(|L|"^)), as L ^ 00 , for j = l,2l 

U{1 — a,l,—djL) = {—djL)~^^~^'>{l + 0{\L\~^)), as L ^ 00 , for j = l,2. J 

Using the identities from (1^51) we obtain 

A/'(a,l,d,L) = £|i^(-d 2 L)-(i+“)(l + 0(|i|-i)), as A ^ 00 , for j = 1,2 1 

U'{1 — a,l,—djL) = (1 — a)(—-I-0(|L|“^)), as L —>■ 00 , for j = l,2. j 
Using (ED: (ED and the relation lim = 0 in the expression (1^^ ^ for ^ 1 (^ 2 ), we obtain 

L—¥oo 

Zi{w 2 ) = o 

which can be written as 21(^2) ~ ^^1^/(0, 1 ,^2), where Ci = 0 {e'^^^(—(see (ED)- 
similar arguments it follows that 


(31) 


(32) 


Using 


Zi{w2) ~ CiM{a,l,d2L) 

Z2{w2) ~ C2M{aA,d2L) 

z[{w 2 ) ~ Ci-M'(a, 1, (i 2 T) = C'iaM(a -I-1,1-I-1, c? 2 U) 

-^ 2 (^ 2 ) ~ C2M'{a,l,d2L) = C2aM{a + 1,1 + ^,d2L) 


for L — 00 , see (l27l) for the dependence of A, di and L of the coefficient Ci, Di, C 2 and -D 2 . Thus, using the 
above asymptotic results for the coefficient Ci, Di, C 2 and D 2 and the asymptotic results for the confluent 


hypergeometric functions given in (1^ and (ED, we get 


Az[{w2) +Bz'^{w2) 

lim —— r- - —- 

L^oo Azi[W2) + BZ2[W2) 


' 1-SfcA 


) lim 

/ r V 


7Vf(l + a, 2, d 2 L)[CiA + C 2 B) 


L^oo M(a, 1, d2L){CiA + C 2 B) 

pi^(-d2T)-(^+“) 

■1^) lim 
- ' L-)-oo 


r(i) 

r(a) 


\-d 2 L)- 


0. 


(34) 
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Substituting this in (0^1) . we obtain 


lim 

L—^co 


/;(0;A) 

/(0;A) 


= k. 


Therefore, equation (E51) becomes fio{0)k = j3\{k)X — I32{k). Using a 
coefficients from ([71)2, we obtain 


t/o/A and expressions for the 


(/i2 - Mo(0)) Tk^ 

(/^2 + Mo( 0)) (/r 2 +/ro(0)). 


which is the classical formula for Saffman-Taylor instability of the leading interface. Similarly, we can also 
recover the the classical formula for Saffman-Taylor instability of the trailing interface by reversing the 
shooting technique as discussed at the end of section 13.11 


Second case (When T —>■ 0): Similar to the previous case, we will first need to get asymptotic approx¬ 
imations for Zi{w 2 ), Z 2 iw 2 ), z'i{w 2 ), and 212 (^ 2 ) in this limit. Notice that in this case, singularities of the 
confluent hypergeometric function of the second kind will arise. Now, we give the following asymptotic 
results from Abramowit z and Stegun [1] 


U{1 - a,l,-djL) = --pYW—r (ln(|djL|)-I-^(1 - a))-I-O(LlnL) 

i (i — aj 

[/(2-a,2,-d,T) = £||5l||d,T|i-2 + 0(lnL), 


(35) 


where we recall that di and d 2 are defined by dsni). Similar to the calculations of the previous case L —>■ 00 , 
we present the dominant terms of the left hand side of (l28l) . It is worth pointing out that due to (l35|) . the 
derivative of the confluent hypergeometric function of the second kind is dominant. From the definition of 
the coefficients Ci, Hi, C 2 and D 2 given in (l27l) and the asymptotic results presented in (IMI) . we obtain 


We remark that 


1 ^1-fcA 

zi{d 2 L) ^ tu{l,2}(, 2 

z[{d2L)r^C, 


r(i) 

r(l-a) 


\diL\-^ 


_ kX \ 

2 2 / 


l-fcAA 

2 ) 


/ 1 _ 1 

\Wl\ Wl\ 

r(l-ha) 


r(i) 

r(2 - a) 


\d2L\-^ 


(36) 


and therefore the asymptotic result for ^((^ 2 ) follows from the definition of the coefficient Ci and Hi, see 
(HZl). From the forms of Ci, Hi, C 2 and H 2 , we get 22 ( 1 ^ 2 ) = 0(21 ( 1 ^ 2 ))• Therefore 


Azi{w2) + Bz2{w2) 
Azi{w2) 


1 . 


Similarly, we obtain 


4(w^2) ~ -H 2 e“=C/'(l - a, 1,- 1 ^ 2 ) 

1 f/(2_a,2,-w;2)e“= 


W{1,2} \ 2 

1 fl-kX' 

IU{1,2} V ^ 

1 _ 1 1 

r(l-a) W{1,2} \d 2 L\ 


r(i) 

r( 2 -a) 


M 2 TI- 


(37) 


( 38 ) 
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Using (1551) . (1571) . and (1551) . it follows that 


L^o /(0;A) 


i-i>0 Azi[W2) + Bz2{'W2j 
Az[(w2) + BZ 2 IW 2 ) 


= k [ 1 — 2 lim 

L^O 


= k \ 1 — 2 lim 

L-s-O 


Azi (W2) 

(fit -m) i+^ (r(i^) (f^ 


= k \ 1 — 2 lim 

L^O 


^(rU^) (f^) 


= A: I 1 — 2 lim 

L-^O 


+b|^' 

V 1*1/ 1*1 


From the definition of di and d 2 given in (l30)) we obtain di/d 2 = fioi—L)/fio{0) and therefore 


lim^ = fc 
L^o /(0;A) 


1 _ 0 ( ( Mo(0) _ r,B fJ-oj-L) 

2 MO) ) ^ MO) 


where A = ^Q{—L)e^^ and B = — ^ '^i(fc)-^+«2 W fc^ol l) ^ ^kL ^ follows that 

lim/ro(0) = [ai(fc) - k'^MiO) - M-L))] X + a2{k). 


L^O 


Using this in equation (1751) . we obtain 

[ai(fc) - fc^(/io(0) - /io(-U))] A + a 2 {k) = l3i{k)X - M^) 

which is equivalent to 


((ai(fc) - mm - k'^i.MO) - M-L))) ^ = -Mk) - Mk)- 

After substituting the values of ai(A:), 02 (^), /3i(fc) and (32{k) and simplifying we obtain 

+ k'^M - ^ 2 ) A = -kM + M- 


Therefore, 


(j = — 


jS + T) , (a^2 - M 

{fJ'2+M ° (//2+Mi)’ 


which is the formula for the growth rate of an interface with surface tension (S' + T), which is what should 
be expected in this limit. Thus we recover the classical formula for the growth rate in this limit. 


4 Conclusions 

We converted a non-standard eigenvalue problem arising in the linear stability analysis of a three-layer Hele- 
Shaw model of enhanced oil recovery to a boundary value problem for Rummer’s equation when the middle 
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layer has a linear viscous profile. We presented the general solution in terms of Frobenius series and discussed 
the convergence properties of these series solutions. We also formally gave the dispersion relation implicitly 
through the existence criterion for non-trivial solutions. In order to recover the well-known physical solutions 
for some limiting cases, we rewrote the general solutions using a different set of fundamental solutions and 
analyzed these for those limiting cases: (i) when the viscous profile of the middle layer approaches a constant 
viscosity, both in the case of a fixed-length middle layer and also as the length of the middle layer appraoches 
infinity; and (ii) when the length of the middle layer approaches zero. We showed that we were thus able to 
recover the correct physical solutions. 
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Appendix: Kummer’s Equation 


Rummer’s equation has the general form 


(fz .dz 

w-—^ -I- (o — w)- - az = U, 

dw-^ dw 


where 6 = 1 and a = ^(1 -I- The two linearly independent solutions are zi{a, 6, w) and Z 2 {a, 6, w) where 
the general expression for 2 : 1 ( 0 , 6, w) is given by 


where 


zi(a,b,w) 


60 


aiw 


9 rj 

a2W CLnW 

622! 633! +■■■+ bnUl 


oo = 1, oi = a. On = a(a + l)(a + 2) ■ ■ ■ (a + n — 1), for n = 2, 3, • • • 

60 = 1, 61 = 6 , bn = 6(6 -|- 1)(6 -I- 2) ■ • • (6 -I- n — 1), for n = 2,3, • • • 

The linearly independent second solution 22 ( 0 , 6 , 10 ) is similarly given by a series which can be easily con¬ 
structed by the method of Frobenius. 
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